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In a recent letter [P.K. Shukla and B. Eliasson, Phys. Rev. Lett. 108, 165007 (2012)] the 
discovery of a new attractive force between protons in a hydrogen plasma was reported that would 
be responsible for the formation of molecules and of a proton lattice. Here we show, based on ab 
initio density functional calculations and general considerations, that these predictions are wrong 
and caused by using linearized quantum hydrodynamics beyond the limits of its applicability. 
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I. INTRODUCTION 

Recently Shukla and Eliasson published in Physical 
Review Letters [l[ a calculation where they applied Lin- 
earized Quantum Hydrodynamics in a simple way to ob- 
tain the effective potential between classical protons in a 
dense hydrogen plasma containing degenerate electrons. 
Their result indicated a novel attractive force between 
two protons. If well-founded, this is not only an interest- 
ing results, but the relatively simple calculation method- 
ology is also potentially one of considerable importance 
as a theorist's tool. It is, therefore, important to discuss 
the implications of these results for real physical systems, 
such as dense quantum plasmas, and to rigorously verify 
the validity and limitations of the method against some 
more fundamental calculation. However, presenting such 
a result cannot be adequately done within the space lim- 
its of a Comment in Physical Review Letters, so it is 
being presented here as a full paper. 

Shukla and Eliasson (SE hereafter) claim [l[ to have 
calculated their potential as that between classical pro- 
tons in a dense hydrogen plasma containing degenerate 
electrons. The hydrogen plasma is considered in equi- 
librium (essentially in the ground state, T = 0), and 
the authors stress that the observed "novel attractive 
force" is different from Friedel oscillations. Furthermore, 
from the existence of a "hard core negative potential", 
Eq. ([1]), they predict the formation of bound states, a 
proton lattice, as well as "critical points and phase tran- 
sitions at nanoscalcs" . If correct, these results could have 
far-reaching consequences for hydrogen and dense quan- 
tum plasmas, in general. 

It is, therefore, of interest to consider the results of 
Shukla and Eliasson more in detail, which is the goal of 
this paper. We first summarize their approach - Lin- 
earized Quantum Hydrodynamics (LQHD). Wen then 
re-evaluate their final formulas transforming to com- 
mon length and energy scales. We then compare their 
screened proton potential to results from ab initio density 
functional (DFT) simulations and observe qualitative de- 
viations. From this, we have to conclude that the LQHD 
results of Ref. [l| are not applicable to dense quantum 



plasmas ad do not provide evidence for attractive inter- 
actions between protons. We conclude by analyzing the 
origin of the deviations of the results of Ref. [l| from 
DFT which are traced to violation of the applicability 
range of linearized quantum hydrodynamics. 

Dense correlated quantum plasmas are presently of 
high interest in many fields including condensed mat- 
ter, astrophysics and laser plasmas, e.g. 0, H|. De- 
spite remarkable theoretical and experimental progress 
over recent decades, even for the simplest plasma sys- 
tem - hydrogen - still interesting questions remain un- 
solved, including details of the phase diagram and the 
behavior under high compression, e.g. 043] ■ The modi- 
fication of the pair interactions between the ions by the 
surrounding plasma is of prime importance for the theo- 
retical understanding of these systems. In contrast to the 
familiar Coulomb potential, 4> l (r) — Q/r, of an ion ob- 
served in vacuo, in a plasma, correlation and quantum ef- 
fects cause screening. At weak non-ideality this gives rise 
to an isotropic Yukawa-type potential, 4>l(r) = ^e~ r /'% 
where l s is the screening length given, in the limit of a 
classical high-temperature plasma, by the Debye radius 
or, in a high-density quantum plasma, by the Thomas- 
Fermi length Ltf- More general screened potentials are 
successfully computed using linear response theory to ob- 
tain the dynamic dielectric function, e'(k, uj) (longitudi- 
nal density response), giving rise to Q 
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where the longitudinal dielectric function derives from 
the dielectric tensor via e'(k, w) = J2ij hkj£ij(k, uj)/k 2 . 
The potential ((T|) typically decays slower than the expo- 
nential Yukawa potential. 

Moreover, screening effects are known to give rise to 
an attractive region of the potential ( "anti-screening" ) in 
nonequilibrium situations such as a charge embedded into 
a streaming flow of oppositely charged particles, which is 
well studied theoretically |9Hll|. and the resulting at- 
tractive force has been confirmed e.g. in dusty plasma 
experiments with streaming ions. No attractive forces 
have been seen in simulations in equilibrium plasmas 
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with two exceptions: i) The formation of bound states 
(such as hydrogen molecules or molecular ions) evidently 
corresponds to a net attractive force between the con- 
stituents (hydrogen atoms, protons), ii) Oscillatory po- 
tentials with (very shallow) negative parts (Friedel os- 
cillations) emerge in a strongly degenerate Fermi gas as 
a consequence of the Fermi edge singularity (essentially 
due to the step character of the zero temperature momen- 
tum distribution of an ideal Fermi gas), e.g. [l2j]. Friedel 
oscillations have been observed in experiments probing 
surface states at very low temperature, e.g. [l3j . 

II. LINEARIZED QUANTUM 
HYDRODYNAMICS (LQHD) 

In Ref. [J SE compute the potential ((XJ) of a proton in 
a hydrogen plasma by 

1. using zero temperature classical hydrodynamic 
equations for the electron density and mean veloc- 
ity, coupled to the Poisson equation for the electro- 
static potential, 

2. adding, in the momentum equation, quantum 
diffraction effects via the Bohm potential Vb (quan- 
tum hydrodynamics, QHD), 

3. adding the pressure of the ideal Fermi gas at zero 
temperature, 

4. adding exchange and correlation effects via an addi- 
tional potential V xc using a simple parametrization, 
Ref. 33 in Q, 

5. neglecting dynamic effects, e'(k, u>) — > e l (k, 0), 

6. solving the resulting hydrodynamic equations in 
linear response for the dielectric function e'(k, 0; r s ) 
(D, in their notation), 

which paramctrically depends on the coupling parame- 
ter (Brueckner parameter) r s = f/as that completely 
defines the plasma state at zero temperature. [Here we 
introduced the mean interparticle distance, f, which is re- 
lated to the unperturbed density by 47rf 3 /3 = n" 1 ]. For 
the explicit form of the used LQHD equations we refer to 
Ref. [l| , for the discussion of V xc and earlier related work, 
see Ref. (TiJ . In the following we restrict ourselves to the 
linearized equations used by SE, i.e. to the LQHD. The 
full QHD is beyond the scope of this paper, for references, 
see Ref. [l|. 

Using LQHD, SE evaluate the static dielectric func- 
tion, e'(k, 0), and screened proton potential, <fi\, at zero 
temperature, as a function of a single parameter a. We 
plot a as a function of the Brueckner parameter in 
Fig. [U SE observe that the screened potential of a pro- 
ton, Eq. (JT|), develops a negative minimum for a > 0.25. 
As can be seen in the figure, a exceeds the value 0.25, in 
a range of coupling parameters 0.61 < r s < 26.22. Let 
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FIG. 1. Coupling parameter a versus Brueckner parameter 
TV An attractive proton potential is predicted in the density 
interval 26.22 > r a > 0.61 at zero temperature. The shaded 
area denotes the range of r s where the plasmon energy is 
smaller than the Fermi energy, &j pG < fceTp (weak coupling). 
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FIG. 2. Screened proton potential of SE for three densities. 
A single negative minimum is observed, shown more clearly 
in the inset. The location of the minimum and its depth are 
shown in Figs. and 01 respectively. Thin vertical lines indi- 
cate the equilibrium nearest neighbor distance of two protons, 
cf. Fig. [3] The black arrow marks the range of values of r 
shown in Fig. [5] 



us now analyze the results for the screened proton poten- 
tial in the whole range of r a values where an attractive 
hydrogen interaction is observed, as done in Ref. [l[. In 
Fig. [5] we plot the potential energy corresponding to the 
screened proton potential in units of eV for three densi- 
ties corresponding to r s = 7,1.5,0.61. At low and high 
densities (i.e. for cases where a < 0.25), we confirm 
that the potential is purely repulsive. In agreement with 
Fig. [U for a > 0.25 a negative minimum develops at a 
distance of several Bohr radii from the proton. The posi- 
tion of the minimum is plotted versus r s in Fig. [3] There 
exists a minimum of this distance of about 5as for r s ~ 1 
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FIG. 3. Position of the attractive potential minimum of SE, 
Ref. PJ, in units of the Bohr radius as a function of r s , 
compared to the mean interparticle distance. Note that the 
ground state nearest neighbor distance of two protons in the 
molecular phase (r s > 2 . . . 3) is in the range of 1.5 . . . 2as- 




FIG. 4. Depth of the attractive potential minimum of SE, 
Ref. [l[ . The deepest value is observed for r s ~ 3 correspond- 
ing to a « 0.6, cf. Fig. Q] 



with a slow (rapid) increase for lower (higher) densities. 
Finally, we analyze the depth of the attractive potential 
minimum. These values are plotted in Fig. |4] for r s be- 
tween 7 and 0.6 below which the minimum vanishes. The 
deepest minimum is observed around r s = 3 and amounts 
to about 6meV corresponding to a temperature of about 
65K. 



III. COMPARISON OF LQHD TO AB INITIO 
DFT RESULTS 



We now compare the LQHD results of SE to what is 
known about dense hydrogen. There have been numer- 
ous studies of hydrogen plasmas based on analytical ap- 
proximations and first-principle simulations as well, e.g. 



& Q- There is hardly any question that the low- 
density ground state consists of -ff 2 - molecules (either in a 
gas, liquid or solid phase). When the density is increased, 
the molecular binding is reduced until, eventually around 
r s ~ 2 ... 3, molecules break up. Even though the de- 
tailed scenario maybe quite complex, e.g. [H, @, involve 
intermediate states (including ff-atoms) and vary with 
temperature it is clear that, around some critical density 
Um (Mott density), all bound states will break up into 
electrons and protons, e.g. [l5j . This pressure ioniza- 
tion is a gradual process and a reasonable estimate for 
Um corresponds to r s w 1.2 ... 1.5. Upon further density 
increase one expects formation of a liquid and, eventu- 
ally, solid phase of protons embedded into a Fermi gas of 
electrons, e.g. [l5|, [l||- While there remain interesting 
questions about the precise values of the different phase 
boundaries [l7j j . the general structure of the hydrogen 
hase diagram at high pressures is well understood, e.g. 



As an illustration, and for direct comparison with SE, 
we present results for the screened proton potential and 
for the interaction of hydrogen atoms embedded in a jel- 
lium background, as obtained from density-functional to- 
tal energy calculations. The ah initio scheme of DFT 
rests on firm theoretical grounds and has been thor- 
oughly tested in recent decades, in particular in appli- 
cation to dense plasmas, see e.g. refs. @,[23. It provides 
a fully consistent treatment of quantum effects and uses 
the fundamental Coulomb interaction between charged 
particles as input. The only uncontrolled approximation 
contained in DFT total-energy calculations is the approx- 
imation applied to the exchange-correlation energy func- 
tional, but the accuracy and limitations of the different 
approximations available have been investigated in detail 
for many applications in solid state physics and molecular 
chemistry [18| ■ Therefore, DFT total-energy calculations 
can serve as a benchmark for the linearized quantum hy- 
drodynamic model of SE. Before proceeding with our own 
calculations we point out that detailed DFT calculations 
have been carried through for H and H2 in jellium be- 
fore, see, e. g., the work by Almbladh et at, Ref. 
by Bonev and Ashcroft, Ref. [2(|, and by Song, Ref. 
Our results concur with their findings. For other recent 
DFT results for dense hydrogen, see Refs. 0,112; HH an d 
references therein. 

In the following we restrict ourselves to spin- 
unpolarized simulations. It is a well-known artefact 
of the usually applied approximations to the exchange- 
correlation energy functional that in vacuo the hydrogen 
molecule undergoes a transition from a spin-unpolarizcd 
to a spin-polarized ground state when the separation of 
the two protons is increased [U,[2^]. However, in their 
spin-polarized DFT-calculations for hydrogen atoms im- 
mersed in jellium, Nazarov et al. [26| have found a spin- 
unpolarized ground-state for 3 < r s < 14. For this rea- 
son, and as no spin-polarization is to be expected at even 
higher electron densities, we have carried through only 
spin-unpolarizcd relaxations for r s — 1.5 and r s = 7. 
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FIG. 5. (color online) Electrostatic potential around an H- 
atom immersed in jellium. r is the distance from the proton. 
DFT data have been calculated in a cubic box (size 15<zo) 
using the PBE-GGA for the exchange-correlation energy- 
functional, a single k-point, and a plane- wave cutoff-energy 
of 200 Ry (r s = 1.5) or 300 Ry (r s = 4 and r a = 7). The 
DFT data (symbols) are compared to the electrostatic poten- 
tial from LQHD of Shukla and Eliasson, Ref. [l[ (lines). The 
black arrow marks the voltage range shown in Fig. 
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FIG. 6. (color online) Interaction energy of two H-atoms 
immersed in jellium for two densities. While for r s = 7 a 
minimum around the hydrogen molecule bond length in vacuo 
(1.4 bohr) is observed, for r s — 1.5, the molecular bound state 
is unstable. The DFT data have been calculated in a cubic 
box (with a size as noted in the inset) using PBE-GGA or 
LDA for the exchange-correlation energy-functional, a single 
or 6 k-points in the irreducible part of the Brillouin zone, and 
a plane- wave cutoff-energy of 100 Ry. 

This is consistent with Ref. f20j . 

The data in Figs. [5] and |5] have been computed with 
the fhimd code [27|. A 1/r Coulomb potential (i.e. no 
pscudopotcntial) has been taken for the proton-electron 
interaction. The cut-off energy of the plane-wave basis- 
set has been chosen at 100 Ry for the calculation of the 



H-H interaction energy, while larger values have been 
used for the calculation of the screened potential of an 
H-atom, or proton, in jellium. The two hydrogen atoms 
forming the H2 molecule are put into a cubic super-cell 
with the length of the edge equal to 25 bohr (15 bohr), in 
case of r s = 7 (r s = 1.5). The super-cell is repeated pe- 
riodically in all three directions. Calculations have been 
carried through using only a single special k-point. Addi- 
tional computations with 6 s pec ial k-points derived from 
a 4 x 4 x 4 back-folding [27], [28[ allow an error estimate, 
as can be seen in Fig. [S] The smearing of the Fermi 
distribution is small, k-g,T = 0.025 eV or 0.1 eV, and 
the total energy has been extrapolated to zero temper- 
ature under the usual assumption of a linear heat ca- 
pacity [27| . The generalized gradient approximation by 
Perdew and Ernzerhof (PBE96) [29| and, for compari- 
son, the local-density approximation, have been applied 
to the exchange-correlation energy functional. 

Let us now consider the results. As a typical low- 
density case, we show in Fig.[5]thc electrostatic potential 
around a proton in jellium for r s = 7. It exhibits a clear 
negative minimum around a distance of r m ,„ w las- 
The minimum becomes shallower with increasing electron 
density, cf. curve for r s =4, and vanishes for r s = 1.5. 
There are further Friedel oscillations at larger distances, 
which, however, are often not well resolved due to com- 
putational limitations and the small size of the super-cell, 
but they are of no relevance for hydrogen at these condi- 
tions, see below. Our data can now be directly compared 
to the LQHD results of SE, shown by lines in Fig. [5] 

In transferring the Shukla-Eliasson results from Fig. [2] 
to Fig. [5] one should note the large changes in scale and 
range. The horizontal axis in Fig. [5] extends out only to 
Qas, while that of Fig. [5] extends more than three times 
further. Also the vertical scale of the potential energy in 
Fig. [2] is roughly six times smaller than that of Fig. [5] 
Thus, on the scale of Fig. (3J the LQHD minima, which 
only appear clearly in the inset of Fig. [3J involve such 
small potential changes on the scale of Fig. [5] that they 
cannot be seen. The small dip in the LQHD potential 
for r s = 1.5 is at about 5Aas near the right hand edge 
of the r-axis, while the dip for r s = 7 is even further 
away, far beyond the right hand edge. In any case, in 
the region of the DFT minimum, the LQHD potential 
is monotonically decreasing, i.e. without minima, and 
far from that of the full calculation, as denoted by the 
symbols. 

In part we ascribe the minimum of the electrostatic 
potential from the DFT calculation in Fig. [5] to the for- 
mation of an H~ ion, i. e. a bound state which Almbladh 
et al. [lj| have found for r s > 1.9. As pointed out by 
Bonev and Ashcroft [2(j, this formation of a negative 
ion sheds doubts on the adequateness of jellium-type (or 
one-component plasma) models for the description of a 
nonidcal hydrogen plasma. The recent ab initio simula- 
tions of hydrogen by Morales et al. Q and Bonev ct al. 
[3] are in fact based on DFT-Born-Oppenheimer molec- 
ular dynamics simulations of ensembles or Monte Carlo 
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simulations and thus fully account for the dynamics of 
the protons. 

The interaction energy of two hydrogen atoms im- 
mersed in a jellium background with r s = 7 or r s = 1.5 is 
displayed in Fig. [6l In fact, the H-H interaction energy 
has been computed by Bonev and Ashcroft [2(| before, 
using the DFT total-energy program VASP. They point 
out that linear response theory would be inadequate for 
this purpose. Furthermore, they describe the fate of the 
hydrogen bond when the electron density is increased: 
For r s > 3 a local energy minimum at a H-H separation 
larger, but still close to, the bond length of the hydrogen 
molecule in vacuum, and an energy barrier towards H2 
association are observed [20|. Our DFT results summa- 
rized in Fig. [6] show the same behavior: For r s = 7 we 
obtain a stable H2-bond at a bond length comparable to 
the bond length of a hydrogen molecule in vacuo, but 
with a much smaller binding energy. This result is not 
sensitive with respect to the approximation applied to 
the exchange-correlation functional. Also the more tech- 
nical effect of the restricted k-point sampling can be read 
from the Figure. Comparing again to the LQHD results 
of Shukla and Eliasson, it is obvious that the interaction 
energy minimum found in the ab initio data is not in- 
cluded in their approach. Thus, the SE screened proton 
potential has, for small densities (e.g. r s — 7), a qualita- 
tively incorrect shape, in particular it completely misses 
the molecular bound state. 

For larger density, i. e. for r s = 1.5, Fig. shows that 
the minimum in our simulations becomes very shallow 
or vanishes at all, within the accuracy of the numeri- 
cal DFT computations. Similarly, we do not resolve ex- 
tremely shallow (below lOOmcV) binding potentials at 
larger distances from the proton, such as the ones re- 
lated to Friedel oscillations. These oscillations could be 
reproduced by DFT using a substantial computational 
effort. However, there is no need for this since such ex- 
tremely small features of the potential are expected to 
be of no relevance for all dense plasma applications since 
they always encounter a finite electron temperature of at 
least one cV. 



IV. DISCUSSION. FAILURE OF LQHD 

The approach by SE in Ref. Q starts from the as- 
sumption that the ions are immobile and embedded into 
a neutralizing background (see Eq. (3) in Ref. pj). We 
have thus used DFT, together with the LDA or the PBE- 
GGA 24[ applied to the exchange-correlation energy 
functional, to calculate the interaction energy of two pro- 
tons immersed in a jellium background, without any fur- 
ther approximation applied to the quantum-mechanical 
problem. For the reliability of the DFT as compared 
to Quantum Monte Carlo simulations we refer e.g. to 
Ref. [23]. Thus the DFT results can serve as a reference 
to evaluate the accuracy of the interacti on p otential as 
derived from the LQHD approach by SE [30j . 



As shown from the comparison to ab initio DFT sim- 
ulations, the SE potential completely misses the bound 
states of protons in low temperature hydrogen. When 
the density is increased the deviations between LQHD 
and DFT potentials become smaller but are still notice- 
able. 

Yet even at higher densities where no molecules exist 
the SE potential is qualitatively wrong because it does 
not show an attractive minimum for r s < 0.61 at all. 
Dielectric theories of the electron gas and electron liquid 
including correlations in the frame of local field correc- 
tions, however, confirm that Friedel oscillations persist 
up to high densities, e.g. [3l|, [HI- The present LQHD 
model of SE also misses the Friedel oscillations, as the 
authors themselves underline. 

Furthermore, we consider two limiting cases of the 
LQHD potential discussed by SE [l[: first, in the limit 
a — > SE recover a Yukawa potential with the screening 
length Ltf- This is correct for high densities, r s — > 0, 
but not for low densities, r s >• 1 where a — > as well, 
cf. Fig. [T] Second, SE recover, for a 3> 1, the exponen- 
tial cosine-screened Coulomb potential [l|. This limit is 
questionable since a cannot exceed 0.65, cf. Fig.[TJ These 
failures rule out any reliable predictions such as novel po- 
tential minima. Besides, the extremely low value of the 
associated binding energy (cf. Fig. 2|) would require a 
particularly accurate theory and careful verification. 

Finally, let us analyze possible reasons for the failure 
of the linearized QHD model in application to dense two- 
component plasmas. 

i. : As any fluid theory, LQHD cannot resolve distances 

below some cut-off r* . This has been pointed out 
by Manfredi and co-workers [l4| who performed a 
test of a ID version of nonlinear QHD for thin metal 
films [3. 

ii. : Shukla and Eliasson note [l| that QHD is valid only 

for Huj pc < fesTp, i.e. for weak coupling. This range 
is indicated by the shaded area in Fig. Q] Since 
the present model of LQHD contains exchange and 
correlation contributions, this range may actually 
be larger [13 ]. 

iii. : Linearized QHD fails when the linearization condi- 

tions are violated, i.e. when the perturbed density 
is comparable with the unperturbed density which 
is common in the description of strong attractive 
interactions. 

While item (ii) is presently open and requires further 
analysis, our DFT simulation results allow us to directly 
verify (L). In a hydrodynamic description of classical 
plasmas the smallest length scale that can be resolved is 
the Debye length. Similarly, in the quantum case, this 
cut-off is expected to be of the order of the Thomas-Fermi 

length, Ltf = frkp / (y/3m e u)p C ) = (n / 12) 1 / 3 rl^ 2 0,3- In 
fact, re-examining the DFT and LQHD data for the 
screened proton potential in Fig. [5] we observe that the 
most dramatic deviations occur on length scales smaller 
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than r* « (3 . . . 4)Ltf- Thus we confirm that LQHD can- 
not, by construction, yield potentials with atomic-scale 
resolution in a dense quantum plasma [33| . 

Yet the most severe limitation is, apparantly, related 
to the linearization of the QHD equations, item (iii). In 
fact, similar linear response calculations for a proton in 
a degenerate electron gas have been done long ago. Alm- 
bladh et al. have performed density functional cal- 
culations and compared them to a linearized version of 
DFT (which, by construction, is essentially more accu- 
rate than any hydrodynamic approach). They observed 
a complete failure of the linear theory in the prediction 
of the electron density screening the proton. Not only is 
the electron density at the proton a factor 2 to 33 (for 
r s = 1 and r s = 6, respectively) too small, but linear 
theory generally predicted Friedel oscillations to occur 
at much too high distances from the proton. Thus, a 
possible explanation for the potential minimum observed 
in the present linearized QHD, cf. Fig. [21 is that it is a 
trace of the Friedel oscillations (the first minimum) that 
is displaced to higher distances from the proton due to 
the linearization, cf. also Fig. [3] It would, therefore, 
be interesting to test this hypothesis by comparing the 
LQHD results of SE to full nonlinear QHD. 



Let us assume for a moment that the attractive poten- 
tial of SE would be a real effect and consider its implica- 
tions for proton crystallization [l[ . To this end, we plot in 
Fig.[3]thc location of the potential minimum and compare 
it to the mean interparticle distance for a given density. 
It is obvious that there is a striking mismatch. Even if 
the system would be at zero temperature protons could 
not occupy the minimum locations simply because their 
density is several orders of magnitude too high (see also 
Fig. [2] where the mean interparticle distance is marked on 
the SE potential), making such a state energetically im- 
possible. We thus have to conclude, based on DFT simu- 
lations and general considerations about the location and 
low depth of the potential minimum, that the predictions 
of SE of novel attractive forces, novel ion lattices, atoms, 
molecules, critical points and phase transitions in dense 
hydrogen are invalid. 
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